In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
library(sp)
library(raster)
library(ggcorrplot)
library(ggh4x)
options(mc.cores = parallel::detectCores())
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")
theme_set(theme_plot())
# Continuous colors
options(ggplot2.continuous.colour = "viridis")
# Discrete colors
scale_colour_discrete <- function(...) {
scale_colour_brewer(palette = "Dark2")
}
scale_fill_discrete <- function(...) {
scale_fill_brewer(palette = "Dark2")
}
cpue_full <- readr::read_csv("data/clean/catch_by_length_q1_q4.csv") %>%
janitor::clean_names() %>%
rename(X = x,
Y = y)
#> Rows: 1759609 Columns: 22
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): Country, haul.id, IDx, ices_rect, id_haul_stomach, species, haul.i...
#> dbl (14): density, year, lat, lon, quarter, Month, sub_div, length_cm, depth...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
# Summaries density by species group an haul
cpue <- cpue_full %>%
mutate(species2 = species,
species2 = ifelse(species == "cod" & length_cm >= 25, "large_cod", species),
species2 = ifelse(species == "cod" & length_cm < 25, "small_cod", species2)) %>%
group_by(haul_id, species2) %>%
summarise(density = sum(density)) %>%
ungroup()
#> mutate: new variable 'species2' (character) with 3 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, species2)
#> summarise: now 27,897 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables
# Trim data with quantiles
ggplot(cpue, aes(density)) +
geom_histogram() +
facet_wrap(~species2, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cpue <- cpue %>%
group_by(species2) %>%
mutate(dens_quants = quantile(density, probs = 0.999)) %>%
filter(density < dens_quants) %>%
ungroup() %>%
dplyr::select(-dens_quants)
#> group_by: one grouping variable (species2)
#> mutate (grouped): new variable 'dens_quants' (double) with 3 unique values and 0% NA
#> filter (grouped): removed 30 rows (<1%), 27,867 rows remaining
#> ungroup: no grouping variables
# Make wide
wcpue <- cpue %>% pivot_wider(names_from = species2, values_from = density)
#> pivot_wider: reorganized (species2, density) into (flounder, large_cod, small_cod) [was 27867x3, now 9373x4]
head(wcpue)
#> # A tibble: 6 × 4
#> haul_id flounder large_cod small_cod
#> <chr> <dbl> <dbl> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 9.41 77.9 0.628
#> 2 1993:1:GFR:SOL:H20:22:32 6.63 5.13 0
#> 3 1993:1:GFR:SOL:H20:23:31 0.953 2.61 0.00459
#> 4 1993:1:GFR:SOL:H20:24:30 1.89 9.71 0
#> 5 1993:1:GFR:SOL:H20:25:2 16.7 400. 4.78
#> 6 1993:1:GFR:SOL:H20:26:3 13.5 179. 2.95
# Left join in trawl information
nrow(wcpue)
#> [1] 9373
hauls <- cpue_full %>% distinct(haul_id, .keep_all = TRUE) %>% dplyr::select(-density)
#> distinct: removed 1,750,236 rows (99%), 9,373 rows remaining
nrow(hauls)
#> [1] 9373
cpue2 <- left_join(wcpue, hauls)
#> Joining, by = "haul_id"
#> left_join: added 20 columns (year, lat, lon, quarter, country, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,373
#> > =======
#> > rows total 9,373
colnames(cpue2)
#> [1] "haul_id" "flounder" "large_cod" "small_cod"
#> [5] "year" "lat" "lon" "quarter"
#> [9] "country" "month" "i_dx" "ices_rect"
#> [13] "sub_div" "length_cm" "id_haul_stomach" "species"
#> [17] "haul_id_size" "substrate" "depth" "temp"
#> [21] "oxy" "sal" "X" "Y"
cpue2 <- cpue2 %>% drop_na(oxy, temp, depth, sal, substrate, flounder, small_cod, large_cod)
#> drop_na: removed 405 rows (4%), 8,968 rows remaining
# Scale variables
cpue2 <- cpue2 %>%
mutate(depth_ct = depth - mean(depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(depth_sq)) / sd(depth_sq),
depth_sc = (depth - mean(depth)) / sd(depth),
temp_ct = temp - mean(temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(temp_sq)) / sd(temp_sq),
temp_sc = (temp - mean(temp)) / sd(temp),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
sal_sc = (sal - mean(sal)) / sd(sal),
fyear = as.factor(year),
fsubstrate = as.factor(substrate))
#> mutate: new variable 'depth_ct' (double) with 149 unique values and 0% NA
#> new variable 'depth_sq' (double) with 149 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 149 unique values and 0% NA
#> new variable 'depth_sc' (double) with 149 unique values and 0% NA
#> new variable 'temp_ct' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sq' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sc' (double) with 8,679 unique values and 0% NA
#> new variable 'oxy_sc' (double) with 8,678 unique values and 0% NA
#> new variable 'sal_sc' (double) with 8,690 unique values and 0% NA
#> new variable 'fyear' (factor) with 28 unique values and 0% NA
#> new variable 'fsubstrate' (factor) with 5 unique values and 0% NA
pred_grid_1 <- read_csv("data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- bind_rows(pred_grid_1, pred_grid_2)
# Scale with respect to data!
pred_grid <- pred_grid %>%
drop_na(substrate) %>%
mutate(X = X / 1000,
Y = Y / 1000,
depth_ct = depth - mean(cpue2$depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
temp_ct = temp - mean(cpue2$temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
oxy_sc = (oxy - mean(cpue2$oxy)) / sd(cpue2$oxy),
sal_sc = (sal - mean(cpue2$sal)) / sd(cpue2$sal),
fyear = as.factor(year),
fsubstrate = as.factor(substrate))
#> drop_na: removed 280 rows (<1%), 592,760 rows remaining
#> mutate: changed 592,760 values (100%) of 'X' (0 new NA)
#> changed 592,760 values (100%) of 'Y' (0 new NA)
#> new variable 'depth_ct' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sq' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sc' (double) with 7,528 unique values and 0% NA
#> new variable 'temp_ct' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sq' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'oxy_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'sal_sc' (double) with 592,760 unique values and 0% NA
#> new variable 'fyear' (factor) with 28 unique values and 0% NA
#> new variable 'fsubstrate' (factor) with 5 unique values and 0% NA
# Split by quarter
pred_grid_q1 <- pred_grid %>% filter(quarter == 1)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
pred_grid_q4 <- pred_grid %>% filter(quarter == 4)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
# Load models
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/spatial_overlap_cache/html")
# hist(cpue2$depth)
# hist(log(cpue2$depth))
#
# cpue_long <- cpue2 %>%
# pivot_longer(c(flounder, small_cod, large_cod)) %>%
# rename(species2 = name, density = value) %>%
# pivot_longer(c(depth, temp, oxy, sal)) %>%
# rename(env_var = name, env_var_value = value)
#
# ggplot(cpue_long, aes(density)) +
# geom_histogram() +
# facet_wrap(~species2, scales = "free", ncol = 1)
#
# ggplot(cpue_long, aes(env_var_value, density)) +
# geom_point(alpha = 0.3) +
# stat_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) +
# stat_smooth(method = "lm", formula = y ~ x + I(x^2), color = "red", se = FALSE) +
# facet_wrap(env_var~species2, scales = "free", ncol = 3) +
# coord_cartesian(expand = 0)
#
# cpue2 <- cpue2 %>% mutate(fle_presence = ifelse(flounder == 0, "N", "Y"),
# l_cod_presence = ifelse(large_cod == 0, "N", "Y"),
# s_cod_presence = ifelse(small_cod == 0, "N", "Y"))
#
# # Biomass density
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = flounder)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = large_cod)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = small_cod)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# # Presence / absence
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = fle_presence)) +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = l_cod_presence)) +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = s_cod_presence)) +
# facet_wrap(~year)
# Q1
spde_q1 <- make_mesh(filter(cpue2, quarter == 1),
xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans",
seed = 42)
#> filter: removed 3,638 rows (41%), 5,330 rows remaining
# All flounder Q4
spde_q4 <- make_mesh(filter(cpue2, quarter == 4),
xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans",
seed = 42)
#> filter: removed 5,330 rows (59%), 3,638 rows remaining
# Small cod model q1 spatial
mcod_s_q1_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining
sanity(mcod_s_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock -0.691928335 1.02212747 -2.69526137
#> 2 fsubstratehard clay -0.132212753 0.63380026 -1.37443843
#> 3 fsubstratehard-bottom complex -0.125952616 0.65479459 -1.40932644
#> 4 fsubstratemud -0.140027403 0.63378872 -1.38223047
#> 5 fsubstratesand -0.294006302 0.63514592 -1.53886944
#> 6 fyear1994 -0.189885126 0.61380222 -1.39291537
#> 7 fyear1995 0.298079563 0.60713833 -0.89188969
#> 8 fyear1996 0.314750404 0.61276989 -0.88625652
#> 9 fyear1997 -0.669099743 0.61556223 -1.87557955
#> 10 fyear1998 -0.040065644 0.60397679 -1.22383840
#> 11 fyear1999 -0.248289902 0.59798496 -1.42031889
#> 12 fyear2000 0.687378415 0.62234748 -0.53240024
#> 13 fyear2001 0.455272793 0.59180912 -0.70465176
#> 14 fyear2002 1.962095060 0.59529336 0.79534152
#> 15 fyear2003 0.320007017 0.61294634 -0.88134573
#> 16 fyear2004 -0.220180305 0.58766498 -1.37198251
#> 17 fyear2005 1.788929063 0.56681351 0.67799500
#> 18 fyear2006 0.259180987 0.58035478 -0.87829347
#> 19 fyear2007 0.957228921 0.56744971 -0.15495207
#> 20 fyear2008 1.371939968 0.56516783 0.26423137
#> 21 fyear2009 0.976393169 0.56340389 -0.12785817
#> 22 fyear2010 1.472598180 0.57147570 0.35252640
#> 23 fyear2011 0.353731080 0.58291380 -0.78875897
#> 24 fyear2012 0.323141901 0.57780572 -0.80933650
#> 25 fyear2013 1.624300129 0.56495747 0.51700384
#> 26 fyear2014 1.075135510 0.57314313 -0.04820438
#> 27 fyear2015 -0.082199342 0.57483185 -1.20884907
#> 28 fyear2016 -0.073825753 0.57662362 -1.20398728
#> 29 fyear2017 -0.127635067 0.58130161 -1.26696528
#> 30 fyear2018 0.629950035 0.57463000 -0.49630407
#> 31 fyear2019 -0.031397966 0.62355142 -1.25353629
#> 32 fyear2020 -1.173909468 0.62458372 -2.39807106
#> 33 depth_sc 1.162865195 0.08949681 0.98745466
#> 34 depth_sq_sc -1.222660927 0.07417559 -1.36804242
#> 35 temp_sc 0.943382195 0.13622800 0.67638023
#> 36 temp_sq_sc -0.472594618 0.07282391 -0.61532686
#> 37 oxy_sc 0.561042569 0.08513249 0.39418595
#> 38 sal_sc -0.001169043 0.10988511 -0.21653990
#> conf.high
#> 1 1.31140470
#> 2 1.11001293
#> 3 1.15742121
#> 4 1.10217566
#> 5 0.95085683
#> 6 1.01314511
#> 7 1.48804882
#> 8 1.51575733
#> 9 0.53738006
#> 10 1.14370711
#> 11 0.92373909
#> 12 1.90715707
#> 13 1.61519734
#> 14 3.12884860
#> 15 1.52135976
#> 16 0.93162190
#> 17 2.89986313
#> 18 1.39665545
#> 19 2.06940992
#> 20 2.47964856
#> 21 2.08064451
#> 22 2.59266996
#> 23 1.49622112
#> 24 1.45562030
#> 25 2.73159642
#> 26 2.19847540
#> 27 1.04445039
#> 28 1.05633577
#> 29 1.01169515
#> 30 1.75620414
#> 31 1.19074036
#> 32 0.05025212
#> 33 1.33827573
#> 34 -1.07727944
#> 35 1.21038416
#> 36 -0.32986238
#> 37 0.72789919
#> 38 0.21420182
# Large cod model q1 spatial
mcod_l_q1_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining
sanity(mcod_l_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 3.35027430 0.73692195 1.90593382 4.79461479
#> 2 fsubstratehard clay 3.86345204 0.47295798 2.93647144 4.79043264
#> 3 fsubstratehard-bottom complex 3.68218901 0.48661209 2.72844685 4.63593118
#> 4 fsubstratemud 3.87755397 0.47296062 2.95056819 4.80453976
#> 5 fsubstratesand 3.89542852 0.47391863 2.96656508 4.82429196
#> 6 fyear1994 0.42180846 0.51298108 -0.58361598 1.42723290
#> 7 fyear1995 0.17395551 0.51194085 -0.82943012 1.17734113
#> 8 fyear1996 1.05226053 0.51303144 0.04673738 2.05778367
#> 9 fyear1997 -0.43001992 0.51539334 -1.44017230 0.58013247
#> 10 fyear1998 -0.58814899 0.51441941 -1.59639250 0.42009452
#> 11 fyear1999 -0.73612259 0.51068879 -1.73705423 0.26480904
#> 12 fyear2000 -0.44805678 0.53791590 -1.50235257 0.60623901
#> 13 fyear2001 -0.05835698 0.51165985 -1.06119185 0.94447790
#> 14 fyear2002 0.32081027 0.51994811 -0.69826929 1.33988984
#> 15 fyear2003 0.29559541 0.52581232 -0.73497780 1.32616861
#> 16 fyear2004 -1.54152654 0.51295582 -2.54690147 -0.53615161
#> 17 fyear2005 -0.23565112 0.49440049 -1.20465827 0.73335602
#> 18 fyear2006 0.42869636 0.49515684 -0.54179321 1.39918593
#> 19 fyear2007 -0.22043688 0.49335112 -1.18738732 0.74651355
#> 20 fyear2008 0.17717789 0.49054227 -0.78426729 1.13862306
#> 21 fyear2009 0.44018278 0.48883021 -0.51790683 1.39827239
#> 22 fyear2010 0.35962700 0.49485452 -0.61027005 1.32952404
#> 23 fyear2011 0.16164975 0.49879137 -0.81596338 1.13926288
#> 24 fyear2012 -0.45179824 0.50054518 -1.43284877 0.52925229
#> 25 fyear2013 -0.06603709 0.49115060 -1.02867458 0.89660039
#> 26 fyear2014 -0.18103119 0.49683537 -1.15481063 0.79274826
#> 27 fyear2015 -0.13023130 0.49483422 -1.10008855 0.83962595
#> 28 fyear2016 -0.20703546 0.49514267 -1.17749726 0.76342635
#> 29 fyear2017 -0.55937561 0.49580712 -1.53113971 0.41238849
#> 30 fyear2018 -0.46754421 0.49737326 -1.44237789 0.50728948
#> 31 fyear2019 -1.04469284 0.53930616 -2.10171350 0.01232781
#> 32 fyear2020 -1.30245419 0.52725884 -2.33586253 -0.26904585
#> 33 depth_sc 0.75692314 0.06450849 0.63048883 0.88335745
#> 34 depth_sq_sc -0.34635027 0.03128711 -0.40767189 -0.28502866
#> 35 temp_sc 0.69470536 0.09776062 0.50309806 0.88631266
#> 36 temp_sq_sc -0.18487760 0.05099597 -0.28482787 -0.08492734
#> 37 oxy_sc 0.32515528 0.06476164 0.19822480 0.45208577
#> 38 sal_sc -0.10165223 0.07887191 -0.25623833 0.05293388
# Flounder model q1 spatial
mfle_q1_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining
sanity(mfle_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 4.217305679 0.50885474 3.2199687197
#> 2 fsubstratehard clay 3.793488470 0.41832849 2.9735796929
#> 3 fsubstratehard-bottom complex 3.709849071 0.42664643 2.8736374301
#> 4 fsubstratemud 3.924364749 0.41826303 3.1045842668
#> 5 fsubstratesand 3.706693351 0.41924775 2.8849828582
#> 6 fyear1994 -0.118952202 0.40958418 -0.9217224479
#> 7 fyear1995 0.043060162 0.40708786 -0.7548173894
#> 8 fyear1996 0.100803326 0.41182762 -0.7063639814
#> 9 fyear1997 -0.595877590 0.41349979 -1.4063222959
#> 10 fyear1998 -0.731133573 0.40940483 -1.5335522945
#> 11 fyear1999 -0.418026699 0.40426484 -1.2103712192
#> 12 fyear2000 -0.008214615 0.42492808 -0.8410583542
#> 13 fyear2001 0.137378708 0.40808148 -0.6624463029
#> 14 fyear2002 0.408813787 0.41429670 -0.4031928310
#> 15 fyear2003 0.246142842 0.41430827 -0.5658864478
#> 16 fyear2004 -0.969159441 0.40269880 -1.7584345819
#> 17 fyear2005 0.025554985 0.39186036 -0.7424771997
#> 18 fyear2006 0.602214237 0.39269289 -0.1674496787
#> 19 fyear2007 0.351034556 0.38921615 -0.4118150731
#> 20 fyear2008 0.301558233 0.39036135 -0.4635359630
#> 21 fyear2009 0.080400536 0.38958103 -0.6831642497
#> 22 fyear2010 0.662656166 0.39188814 -0.1054304685
#> 23 fyear2011 0.247734357 0.39512575 -0.5266978816
#> 24 fyear2012 0.100854403 0.39259612 -0.6686198511
#> 25 fyear2013 0.466108433 0.38842681 -0.2951941212
#> 26 fyear2014 0.156389468 0.39370883 -0.6152656656
#> 27 fyear2015 -0.267006737 0.39316519 -1.0375963409
#> 28 fyear2016 -0.151282214 0.39280337 -0.9211626677
#> 29 fyear2017 -0.115556315 0.39224170 -0.8843359205
#> 30 fyear2018 0.238284238 0.39501972 -0.5359401781
#> 31 fyear2019 -0.186458600 0.42712577 -1.0236097344
#> 32 fyear2020 -0.327708208 0.41651490 -1.1440624056
#> 33 depth_sc 0.097959908 0.04974479 0.0004619065
#> 34 depth_sq_sc -0.055185737 0.02318019 -0.1006180743
#> 35 temp_sc 0.837214787 0.08491743 0.6707796801
#> 36 temp_sq_sc -0.009634014 0.04451165 -0.0968752532
#> 37 oxy_sc 0.563365911 0.05032856 0.4647237384
#> 38 sal_sc 0.357375668 0.07569706 0.2090121600
#> conf.high
#> 1 5.214642638
#> 2 4.613397247
#> 3 4.546060712
#> 4 4.744145232
#> 5 4.528403844
#> 6 0.683818045
#> 7 0.840937714
#> 8 0.907970633
#> 9 0.214567115
#> 10 0.071285149
#> 11 0.374317821
#> 12 0.824629125
#> 13 0.937203718
#> 14 1.220820404
#> 15 1.058172133
#> 16 -0.179884301
#> 17 0.793587169
#> 18 1.371878154
#> 19 1.113884186
#> 20 1.066652430
#> 21 0.843965321
#> 22 1.430742801
#> 23 1.022166596
#> 24 0.870328657
#> 25 1.227410987
#> 26 0.928044602
#> 27 0.503582867
#> 28 0.618598239
#> 29 0.653223290
#> 30 1.012508654
#> 31 0.650692535
#> 32 0.488645990
#> 33 0.195457910
#> 34 -0.009753399
#> 35 1.003649893
#> 36 0.077607226
#> 37 0.662008083
#> 38 0.505739176
# Small cod model q4 spatial
mcod_s_q4_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining
sanity(mcod_s_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 0.09435785 0.74557893 -1.36694999 1.55566570
#> 2 fsubstratehard clay -0.78082520 0.55629295 -1.87113935 0.30948894
#> 3 fsubstratehard-bottom complex -1.30971300 0.57320481 -2.43317377 -0.18625222
#> 4 fsubstratemud -1.01679675 0.55400257 -2.10262183 0.06902833
#> 5 fsubstratesand -0.81740941 0.55604348 -1.90723460 0.27241578
#> 6 fyear1994 0.09575997 0.69520069 -1.26680834 1.45832828
#> 7 fyear1995 1.08520169 0.65733901 -0.20315911 2.37356248
#> 8 fyear1996 -0.65602062 0.69879535 -2.02563434 0.71359309
#> 9 fyear1997 0.70967734 0.63506435 -0.53502593 1.95438060
#> 10 fyear1998 0.42678207 0.65481121 -0.85662432 1.71018847
#> 11 fyear1999 0.70832897 0.63583473 -0.53788420 1.95454215
#> 12 fyear2000 1.15887477 0.61295487 -0.04249471 2.36024425
#> 13 fyear2001 1.94399364 0.60273914 0.76264662 3.12534065
#> 14 fyear2002 1.48091157 0.60100437 0.30296466 2.65885849
#> 15 fyear2003 -0.17510629 0.61432459 -1.37916036 1.02894778
#> 16 fyear2004 2.43321357 0.59484838 1.26733217 3.59909496
#> 17 fyear2005 1.51499545 0.57814968 0.38184290 2.64814800
#> 18 fyear2006 1.91959803 0.57576174 0.79112576 3.04807031
#> 19 fyear2007 2.22121501 0.57216227 1.09979756 3.34263245
#> 20 fyear2008 1.67400570 0.57090638 0.55504976 2.79296165
#> 21 fyear2009 2.37379404 0.56946394 1.25766523 3.48992284
#> 22 fyear2010 1.72512674 0.57286923 0.60232369 2.84792979
#> 23 fyear2011 0.87635038 0.57624420 -0.25306751 2.00576826
#> 24 fyear2012 1.31752499 0.58228856 0.17626038 2.45878960
#> 25 fyear2013 1.98309676 0.57388746 0.85829801 3.10789550
#> 26 fyear2014 0.88658990 0.58050683 -0.25118259 2.02436238
#> 27 fyear2015 0.50716436 0.58397545 -0.63740649 1.65173521
#> 28 fyear2016 0.16799697 0.58158813 -0.97189482 1.30788876
#> 29 fyear2017 1.49569425 0.57485496 0.36899923 2.62238927
#> 30 fyear2018 0.53595781 0.57927770 -0.59940563 1.67132124
#> 31 fyear2019 0.05818354 0.64153030 -1.19919274 1.31555983
#> 32 fyear2020 0.26687992 0.61174312 -0.93211457 1.46587441
#> 33 depth_sc 0.17234560 0.09341009 -0.01073482 0.35542602
#> 34 depth_sq_sc -1.27936173 0.08813100 -1.45209531 -1.10662815
#> 35 temp_sc 0.76108938 0.13370443 0.49903352 1.02314524
#> 36 temp_sq_sc -0.20875558 0.07211844 -0.35010512 -0.06740604
#> 37 oxy_sc 1.07125991 0.10102534 0.87325388 1.26926595
#> 38 sal_sc 0.44360031 0.11241137 0.22327807 0.66392255
# Large cod model q4 spatial
mcod_l_q4_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining
sanity(mcod_l_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 4.21170576 0.68969740 2.85992369 5.56348782
#> 2 fsubstratehard clay 3.86010397 0.51098291 2.85859587 4.86161206
#> 3 fsubstratehard-bottom complex 3.38036271 0.52245906 2.35636178 4.40436364
#> 4 fsubstratemud 3.77126497 0.50969618 2.77227881 4.77025113
#> 5 fsubstratesand 3.71416028 0.51148436 2.71166936 4.71665120
#> 6 fyear1994 0.42759184 0.56627641 -0.68228954 1.53747321
#> 7 fyear1995 1.17104458 0.56088597 0.07172829 2.27036088
#> 8 fyear1996 0.30666879 0.57086084 -0.81219790 1.42553548
#> 9 fyear1997 -0.27644075 0.55060815 -1.35561290 0.80273140
#> 10 fyear1998 -0.12468170 0.55891250 -1.22013007 0.97076667
#> 11 fyear1999 -0.50545258 0.55787366 -1.59886486 0.58795970
#> 12 fyear2000 0.16521615 0.53777897 -0.88881126 1.21924356
#> 13 fyear2001 0.17695597 0.53395591 -0.86957838 1.22349032
#> 14 fyear2002 0.82122342 0.52673599 -0.21116016 1.85360699
#> 15 fyear2003 -1.41070205 0.54015717 -2.46939066 -0.35201344
#> 16 fyear2004 0.11698731 0.53032500 -0.92243058 1.15640520
#> 17 fyear2005 0.60911272 0.51052745 -0.39150268 1.60972813
#> 18 fyear2006 0.42784518 0.51125268 -0.57419165 1.42988202
#> 19 fyear2007 0.57901828 0.50949887 -0.41958116 1.57761772
#> 20 fyear2008 0.71078769 0.50633284 -0.28160644 1.70318182
#> 21 fyear2009 0.48858518 0.50727701 -0.50565948 1.48282984
#> 22 fyear2010 0.74756835 0.50853799 -0.24914779 1.74428449
#> 23 fyear2011 0.01460874 0.51048187 -0.98591734 1.01513482
#> 24 fyear2012 -0.40274524 0.51964102 -1.42122292 0.61573244
#> 25 fyear2013 0.06254185 0.51097323 -0.93894727 1.06403097
#> 26 fyear2014 -0.13121984 0.51357630 -1.13781090 0.87537122
#> 27 fyear2015 -0.02210225 0.51330778 -1.02816700 0.98396251
#> 28 fyear2016 -0.71613915 0.51517454 -1.72586269 0.29358439
#> 29 fyear2017 -0.41320127 0.51291329 -1.41849284 0.59209030
#> 30 fyear2018 -1.08808137 0.51586811 -2.09916429 -0.07699845
#> 31 fyear2019 -1.21055219 0.56145868 -2.31099099 -0.11011339
#> 32 fyear2020 -1.50211588 0.54522694 -2.57074105 -0.43349071
#> 33 depth_sc 0.37248818 0.07157203 0.23220958 0.51276677
#> 34 depth_sq_sc -0.57573368 0.05693097 -0.68731633 -0.46415102
#> 35 temp_sc 0.28283821 0.09928947 0.08823443 0.47744200
#> 36 temp_sq_sc -0.08858374 0.05410450 -0.19462660 0.01745913
#> 37 oxy_sc 0.73749983 0.07553416 0.58945561 0.88554405
#> 38 sal_sc 0.24195972 0.08730570 0.07084370 0.41307574
# Flounder model q4 spatial
mfle_q4_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining
sanity(mfle_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 2.51946863 0.55705485 1.42766119 3.61127606
#> 2 fsubstratehard clay 2.15626423 0.51660542 1.14373621 3.16879225
#> 3 fsubstratehard-bottom complex 2.00553041 0.52214542 0.98214418 3.02891664
#> 4 fsubstratemud 2.09766143 0.51521436 1.08785984 3.10746303
#> 5 fsubstratesand 2.17656896 0.51672767 1.16380134 3.18933657
#> 6 fyear1994 -0.47801146 0.56314064 -1.58174682 0.62572391
#> 7 fyear1995 0.70747231 0.54561884 -0.36192096 1.77686558
#> 8 fyear1996 0.21787202 0.56110562 -0.88187479 1.31761884
#> 9 fyear1997 0.21476366 0.52933021 -0.82270450 1.25223181
#> 10 fyear1998 -0.27099120 0.53917133 -1.32774758 0.78576518
#> 11 fyear1999 -1.03874014 0.54195660 -2.10095556 0.02347528
#> 12 fyear2000 0.02358672 0.52036980 -0.99631935 1.04349279
#> 13 fyear2001 -0.13931059 0.51746380 -1.15352100 0.87489981
#> 14 fyear2002 0.79548001 0.50897871 -0.20209993 1.79305996
#> 15 fyear2003 -1.03066836 0.50997714 -2.03020519 -0.03113152
#> 16 fyear2004 0.39486395 0.50887822 -0.60251904 1.39224694
#> 17 fyear2005 0.43871414 0.49061216 -0.52286802 1.40029631
#> 18 fyear2006 0.46242684 0.48906494 -0.49612283 1.42097650
#> 19 fyear2007 0.59181644 0.48775073 -0.36415742 1.54779031
#> 20 fyear2008 0.80776983 0.48502260 -0.14285700 1.75839667
#> 21 fyear2009 0.71303600 0.48347509 -0.23455776 1.66062976
#> 22 fyear2010 1.20423026 0.48405021 0.25550928 2.15295124
#> 23 fyear2011 0.75829681 0.48413917 -0.19059854 1.70719215
#> 24 fyear2012 0.36522721 0.48772824 -0.59070257 1.32115699
#> 25 fyear2013 0.68792578 0.48757151 -0.26769682 1.64354837
#> 26 fyear2014 0.34406155 0.48797468 -0.61235125 1.30047435
#> 27 fyear2015 0.37608552 0.48740447 -0.57920969 1.33138074
#> 28 fyear2016 0.25798676 0.48543592 -0.69345017 1.20942369
#> 29 fyear2017 0.47267936 0.48773506 -0.48326379 1.42862251
#> 30 fyear2018 0.25476915 0.48878170 -0.70322538 1.21276368
#> 31 fyear2019 0.26560882 0.52792617 -0.76910747 1.30032511
#> 32 fyear2020 -0.12599992 0.51004667 -1.12567302 0.87367319
#> 33 depth_sc -0.18628399 0.06551021 -0.31468164 -0.05788635
#> 34 depth_sq_sc -0.70612932 0.06538288 -0.83427741 -0.57798124
#> 35 temp_sc 0.64874416 0.09993836 0.45286857 0.84461975
#> 36 temp_sq_sc -0.06611203 0.05554101 -0.17497041 0.04274636
#> 37 oxy_sc 1.20548898 0.07915864 1.05034089 1.36063707
#> 38 sal_sc 0.22610329 0.08758623 0.05443743 0.39776914
# Q1
# small cod
small_cod_q1 <- tidy(mcod_s_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Small cod",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# large cod
large_cod_q1 <- tidy(mcod_l_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Large cod",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# flounder
flounder_q1 <- tidy(mfle_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Flounder",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
# small cod
small_cod_q4 <- tidy(mcod_s_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Small cod",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# large cod
large_cod_q4 <- tidy(mcod_l_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Large cod",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# flounder
flounder_q4 <- tidy(mfle_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Flounder",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
coefs <- bind_rows(small_cod_q1,
large_cod_q1,
flounder_q1,
small_cod_q4,
large_cod_q4,
flounder_q4) %>%
mutate(param_group = ifelse(term %in% c("fsubstratebedrock",
"fsubstratehard clay",
"fsubstratehard-bottom complex",
"fsubstratemud",
"fsubstratesand"),
"Substrate",
"Continious")) %>%
mutate(term = ifelse(term == "fsubstratebedrock", "Bedrock", term),
term = ifelse(term == "fsubstratehard clay", "Hard clay", term),
term = ifelse(term == "fsubstratehard-bottom complex", "Hard-bottom\ncomplex", term),
term = ifelse(term == "fsubstratemud", "Mud", term),
term = ifelse(term == "fsubstratesand", "Sand", term),
term = ifelse(term == "depth_sc", "Depth", term),
term = ifelse(term == "depth_sq_sc", "Depth^2", term),
term = ifelse(term == "temp_sc", "Temperature", term),
term = ifelse(term == "temp_sq_sc", "Temperature^2", term),
term = ifelse(term == "oxy_sc", "Oxygen", term),
term = ifelse(term == "sal_sc", "Salinity", term))
#> mutate: new variable 'param_group' (character) with 2 unique values and 0% NA
#> mutate: changed 66 values (29%) of 'term' (0 new NA)
coefs %>%
filter(!grepl('year', term)) %>%
ggplot(aes(term, estimate, color = species, shape = factor(quarter))) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
labs(shape = "Quarter") +
facet_wrap(~param_group, scales = "free") +
labs(x = "Predictor", y = "Standardized coefficient") +
theme_plot() +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
aspect.ratio = 1) +
NULL
#> filter: removed 162 rows (71%), 66 rows remaining
ggsave("figures/habitat_coefs.pdf", width = 20, height = 20, units = "cm")
mcod_s_q1_pred <- predict(mcod_s_q1_space, newdata = pred_grid_q1)
mcod_s_q4_pred <- predict(mcod_s_q4_space, newdata = pred_grid_q4)
mcod_l_q1_pred <- predict(mcod_l_q1_space, newdata = pred_grid_q1)
mcod_l_q4_pred <- predict(mcod_l_q4_space, newdata = pred_grid_q4)
mfle_q1_pred <- predict(mfle_q1_space, newdata = pred_grid_q1)
mfle_q4_pred <- predict(mfle_q4_space, newdata = pred_grid_q4)
Plot predictions on map
p1 <- plot_map +
geom_raster(data = bind_rows(mcod_s_q1_pred, mcod_s_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Cod < 25 cm") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
p2 <- plot_map +
geom_raster(data = bind_rows(mcod_l_q1_pred, mcod_l_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Cod > 25 cm") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
p3 <- plot_map +
geom_raster(data = bind_rows(mfle_q1_pred, mfle_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Flounder") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
((p1 + p2) / (p3 + plot_spacer())) + plot_annotation(tag_levels = "A") &
theme(legend.text = element_text(angle = 90))
# Check flounder...
plot_map +
geom_raster(data = mfle_q1_pred,
aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt")
plot_map +
geom_raster(data = mfle_q4_pred,
aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt")
all_pred_q1 <- mcod_s_q1_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA
all_pred_q1 <- all_pred_q1 %>%
mutate(large_cod = exp(mcod_l_q1_pred$est),
fle = exp(mfle_q1_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#> new variable 'fle' (double) with 296,380 unique values and 0% NA
all_pred_q4 <- mcod_s_q4_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA
all_pred_q4 <- all_pred_q4 %>%
mutate(large_cod = exp(mcod_l_q4_pred$est),
fle = exp(mfle_q4_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#> new variable 'fle' (double) with 296,380 unique values and 0% NA
all_pred <- bind_rows(all_pred_q1, all_pred_q4)
# Calculate biomass overlap
loc_collocspfn <- function(prey, pred) {
p_prey <- prey/sum(prey, na.rm = T)
p_pred <- pred/sum(pred, na.rm = T)
(p_prey*p_pred)/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}
loc_collocspfn_tot <- function(prey, pred) {
p_prey <- prey/sum(prey, na.rm = T)
p_pred <- pred/sum(pred, na.rm = T)
sum((p_prey*p_pred))/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}
all_pred <- all_pred %>%
group_by(year, quarter) %>%
mutate(scod_fle_ovr = loc_collocspfn(pred = small_cod, prey = fle),
scod_fle_ovr_tot = loc_collocspfn_tot(pred = small_cod, prey = fle),
lcod_fle_ovr = loc_collocspfn(pred = large_cod, prey = fle),
lcod_fle_ovr_tot = loc_collocspfn_tot(pred = large_cod, prey = fle)) %>%
ungroup()
#> group_by: 2 grouping variables (year, quarter)
#> mutate (grouped): new variable 'scod_fle_ovr' (double) with 592,760 unique values and 0% NA
#> new variable 'scod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> new variable 'lcod_fle_ovr' (double) with 592,760 unique values and 0% NA
#> new variable 'lcod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> ungroup: no grouping variables
Plot
plot_map_fc +
geom_raster(data = all_pred %>% filter(quarter == 1),
aes(x = X, y = Y, fill = scod_fle_ovr)) +
scale_fill_viridis_c(trans = "sqrt", name = "small cod-flounder") +
facet_wrap(~ year, ncol = 5) +
theme(legend.text = element_text(angle = 90)) +
NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).
plot_map_fc +
geom_raster(data = all_pred %>% filter(quarter == 4),
aes(x = X, y = Y, fill = lcod_fle_ovr)) +
scale_fill_viridis_c(trans = "sqrt", name = "large cod-flounder") +
facet_wrap(~ year, ncol = 5) +
theme(legend.text = element_text(angle = 90)) +
NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).
# In space for selected year
all_pred2 <- all_pred %>%
filter(year %in% c(1995, 2015)) %>%
dplyr::select(X, Y, scod_fle_ovr, lcod_fle_ovr, quarter, year) %>%
pivot_longer(c(scod_fle_ovr, lcod_fle_ovr)) %>%
rename(overlap = name)
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
#> pivot_longer: reorganized (scod_fle_ovr, lcod_fle_ovr) into (name, value) [was 42340x6, now 84680x6]
#> rename: renamed one variable (overlap)
p1 <- plot_map_fc +
geom_raster(data = filter(all_pred2, overlap == "scod_fle_ovr"),
aes(x = X*1000, y = Y*1000, fill = value)) +
scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
facet_grid(year ~ quarter) +
theme(legend.text = element_text(angle = 90)) +
ggtitle("Small cod : Flonder") +
NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining
p2 <- plot_map_fc +
geom_raster(data = filter(all_pred2, overlap == "lcod_fle_ovr"),
aes(x = X*1000, y = Y*1000, fill = value)) +
scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
facet_grid(year ~ quarter) +
theme(legend.text = element_text(angle = 90)) +
ggtitle("Large cod : Flonder") +
NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining
p1 + p2
ggsave("figures/spatial_overlap.pdf", width = 22, height = 14, units = "cm")
# Over time
p5 <- all_pred %>%
distinct(year, quarter, .keep_all = TRUE) %>%
ggplot(aes(year, lcod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) +
labs(color = "Quarter", x = "Year", y = "Overlap index") +
guides(fill = "none") +
stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
ggtitle("Large cod : Flounder") +
geom_point()
#> distinct: removed 592,704 rows (>99%), 56 rows remaining
p6 <- all_pred %>%
distinct(year, quarter, .keep_all = TRUE) %>%
ggplot(aes(year, scod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) +
labs(color = "Quarter", x = "Year", y = "Overlap index") +
guides(fill = "none") +
stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
ggtitle("Small cod : Flounder") +
geom_point()
#> distinct: removed 592,704 rows (>99%), 56 rows remaining
(p5 | p6) + plot_annotation(tag_levels = "A") & theme(aspect.ratio = 1)
ggsave("figures/spatial_overlap_time.pdf", width = 20, height = 20, units = "cm")
# Scale with respect to data!
nd_depth <- data.frame(depth = seq(min(cpue2$depth), max(cpue2$depth), length.out = 100)) %>%
mutate(depth_ct = depth - mean(cpue2$depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
temp_sq_sc = 0,
temp_sc = 0,
oxy_sc = 0,
sal_sc = 0,
year = 1999L,
fyear = as.factor(1999),
fsubstrate = "mud")
#> mutate: new variable 'depth_ct' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 100 unique values and 0% NA
#> new variable 'depth_sc' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with one unique value and 0% NA
#> new variable 'temp_sc' (double) with one unique value and 0% NA
#> new variable 'oxy_sc' (double) with one unique value and 0% NA
#> new variable 'sal_sc' (double) with one unique value and 0% NA
#> new variable 'year' (integer) with one unique value and 0% NA
#> new variable 'fyear' (factor) with one unique value and 0% NA
#> new variable 'fsubstrate' (character) with one unique value and 0% NA
# Q1
margin_small_cod_pred_q1_depth <- predict(mcod_s_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q1_depth <- predict(mcod_l_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q1_depth <- predict(mfle_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
margin_small_cod_pred_q4_depth <- predict(mcod_s_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q4_depth <- predict(mcod_l_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q4_depth <- predict(mfle_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margins_depth <- bind_rows(margin_small_cod_pred_q1_depth,
margin_large_cod_pred_q1_depth,
margin_fle_pred_q1_depth,
margin_small_cod_pred_q4_depth,
margin_large_cod_pred_q4_depth,
margin_fle_pred_q4_depth) %>%
mutate(species = ifelse(species == "flounder", "Flounder", species),
species = ifelse(species == "small_cod", "Small cod", species),
species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)
# Temperature
nd_temp <- data.frame(temp = seq(min(cpue2$temp), max(cpue2$temp), length.out = 100)) %>%
mutate(temp_ct = temp - mean(cpue2$temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
depth_sq_sc = 0,
depth_sc = 0,
oxy_sc = 0,
sal_sc = 0,
year = 1999L,
fyear = as.factor(1999),
fsubstrate = "mud")
#> mutate: new variable 'temp_ct' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 100 unique values and 0% NA
#> new variable 'temp_sc' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with one unique value and 0% NA
#> new variable 'depth_sc' (double) with one unique value and 0% NA
#> new variable 'oxy_sc' (double) with one unique value and 0% NA
#> new variable 'sal_sc' (double) with one unique value and 0% NA
#> new variable 'year' (integer) with one unique value and 0% NA
#> new variable 'fyear' (factor) with one unique value and 0% NA
#> new variable 'fsubstrate' (character) with one unique value and 0% NA
# Q1
margin_small_cod_pred_q1_temp <- predict(mcod_s_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q1_temp <- predict(mcod_l_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q1_temp <- predict(mfle_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
margin_small_cod_pred_q4_temp <- predict(mcod_s_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q4_temp <- predict(mcod_l_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q4_temp <- predict(mfle_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margins_temp <- bind_rows(margin_small_cod_pred_q1_temp,
margin_large_cod_pred_q1_temp,
margin_fle_pred_q1_temp,
margin_small_cod_pred_q4_temp,
margin_large_cod_pred_q4_temp,
margin_fle_pred_q4_temp) %>%
mutate(species = ifelse(species == "flounder", "Flounder", species),
species = ifelse(species == "small_cod", "Small cod", species),
species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)
# Plot!
p1 <- ggplot(margins_depth,
aes(depth, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
color = species, fill = species)) +
geom_line() +
facet_wrap(~quarter, scales = "free") +
geom_ribbon(alpha = 0.4, color = NA) +
coord_cartesian(xlim = c(10, 170)) +
theme(aspect.ratio = 1) +
labs(x = "Depth (m)", y = "Biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
p2 <- ggplot(margins_temp,
aes(temp, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
color = species, fill = species)) +
geom_line() +
facet_wrap(~quarter, scales = "free") +
geom_ribbon(alpha = 0.4, color = NA) +
coord_cartesian(xlim = c(1, 14)) +
theme(aspect.ratio = 1) +
labs(x = "Temperature (°C)", y = "Biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
(p1 / p2) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")
ggsave("figures/supp/conditional_effects.pdf", width = 20, height = 20, units = "cm")
# Scaled
p3 <- margins_temp %>%
ungroup() %>%
group_by(species, quarter) %>%
mutate(max_est = max(exp(est)),
est_sc = exp(est) / max_est) %>%
ggplot(aes(temp, est_sc, color = species, fill = species)) +
geom_line(size = 1.3) +
facet_wrap(~quarter, scales = "free") +
coord_cartesian(xlim = c(1, 14)) +
theme(aspect.ratio = 1) +
labs(x = "Temperature (°C)", y = "Scaled biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#> new variable 'est_sc' (double) with 595 unique values and 0% NA
p4 <- margins_depth %>%
ungroup() %>%
group_by(species, quarter) %>%
mutate(max_est = max(exp(est)),
est_sc = exp(est) / max_est) %>%
ggplot(aes(depth, est_sc, color = species, fill = species)) +
geom_line(size = 1.3) +
facet_wrap(~quarter, scales = "free") +
coord_cartesian(xlim = c(10, 170)) +
theme(aspect.ratio = 1) +
labs(x = "Depth (m)", y = "Scaled biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#> new variable 'est_sc' (double) with 595 unique values and 0% NA
(p3 / p4) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")
ggsave("figures/conditional_effects.pdf", width = 20, height = 20, units = "cm")
head(all_pred) %>% as.data.frame()
#> X Y depth year lon lat substrate quarter oxy temp
#> 1 450 5984 11 1993 14.23718 54.00188 sand 1 8.713132 2.821855
#> 2 454 5984 11 1993 14.29820 54.00225 sand 1 8.743857 2.811376
#> 3 458 5984 11 1993 14.35922 54.00259 sand 1 8.745697 2.806495
#> 4 462 5984 12 1993 14.42024 54.00290 sand 1 8.759304 2.802808
#> 5 466 5984 12 1993 14.48127 54.00318 sand 1 8.765212 2.805644
#> 6 470 5984 11 1993 14.54229 54.00343 sand 1 8.805174 2.792613
#> sal sub_div ices_rect density_saduria biomass_spr biomass_her
#> 1 7.293935 24 37G4 0 NA NA
#> 2 7.346682 24 37G4 0 NA NA
#> 3 7.383028 24 37G4 0 NA NA
#> 4 7.418652 24 37G4 0 NA NA
#> 5 7.501428 24 37G4 0 NA NA
#> 6 7.585449 24 37G4 0 NA NA
#> biomass_spr_sd biomass_her_sd depth_ct depth_sq depth_sq_sc depth_sc
#> 1 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 2 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 3 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 4 NA NA -35.63035 1269.522 0.4131041 -1.306567
#> 5 NA NA -35.63035 1269.522 0.4131041 -1.306567
#> 6 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> temp_ct temp_sq temp_sq_sc temp_sc oxy_sc sal_sc fyear fsubstrate
#> 1 -3.484890 12.14446 0.5884181 -1.305532 1.178372 -1.098960 1993 sand
#> 2 -3.495369 12.21760 0.5969916 -1.309458 1.190285 -1.082916 1993 sand
#> 3 -3.500250 12.25175 0.6009941 -1.311287 1.190999 -1.071861 1993 sand
#> 4 -3.503937 12.27757 0.6040213 -1.312668 1.196275 -1.061026 1993 sand
#> 5 -3.501101 12.25771 0.6016927 -1.311606 1.198566 -1.035849 1993 sand
#> 6 -3.514132 12.34912 0.6124078 -1.316487 1.214061 -1.010293 1993 sand
#> est est_non_rf est_rf omega_s epsilon_st small_cod large_cod
#> 1 -3.656425 -3.277790 -0.3786354 -0.5617690 0.1831337 0.02582466 5.541801
#> 2 -3.865310 -3.278880 -0.5864303 -0.7847163 0.1982860 0.02095642 5.108481
#> 3 -4.067168 -3.282109 -0.7850593 -0.9984999 0.2134406 0.01712582 4.721942
#> 4 -4.041875 -3.169858 -0.8720176 -1.1006398 0.2286222 0.01756450 4.789478
#> 5 -4.125475 -3.166499 -0.9589759 -1.2027798 0.2438039 0.01615582 4.634219
#> 6 -4.325476 -3.279542 -1.0459342 -1.3049197 0.2589855 0.01322725 4.266585
#> fle scod_fle_ovr scod_fle_ovr_tot lcod_fle_ovr lcod_fle_ovr_tot
#> 1 14.79964 4.175460e-08 0.3650466 9.478820e-07 0.5484333
#> 2 14.55492 3.332310e-08 0.3650466 8.593177e-07 0.5484333
#> 3 14.27225 2.670313e-08 0.3650466 7.788702e-07 0.5484333
#> 4 14.71945 2.824528e-08 0.3650466 8.147639e-07 0.5484333
#> 5 15.16271 2.676235e-08 0.3650466 8.120923e-07 0.5484333
#> 6 15.55304 2.247519e-08 0.3650466 7.669161e-07 0.5484333
all_pred_sub <- all_pred %>%
dplyr::select(small_cod, large_cod, fle, year, temp, oxy, sal, depth, quarter) %>%
pivot_longer(c(small_cod, large_cod, fle)) %>%
rename(Species = name, density = value)
#> pivot_longer: reorganized (small_cod, large_cod, fle) into (name, value) [was 592760x9, now 1778280x8]
#> rename: renamed 2 variables (Species, density)
wm <- all_pred_sub %>%
pivot_longer(c(temp, oxy, sal, depth)) %>%
mutate(name = ifelse(name == "temp", "Temperature (°C)", name),
name = ifelse(name == "sal", "Salinity", name),
name = ifelse(name == "depth", "Depth", name),
name = ifelse(name == "oxy", "Oxygen (ml/L)", name)) %>%
group_by(year, quarter, Species, name) %>%
summarise("Weighted median" = hutils::weighted_quantile(v = value, w = density, p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.75))) %>%
rename(env_var = name) %>%
mutate(Species = ifelse(Species == "fle", "Flounder", Species),
Species = ifelse(Species == "large_cod", "Large cod", Species),
Species = ifelse(Species == "small_cod", "Small cod", Species))
#> pivot_longer: reorganized (temp, oxy, sal, depth) into (name, value) [was 1778280x8, now 7113120x6]
#> mutate: changed 7,113,120 values (100%) of 'name' (0 new NA)
#> group_by: 4 grouping variables (year, quarter, Species, name)
#> summarise: now 672 rows and 7 columns, 3 group variables remaining (year, quarter, Species)
#> rename: renamed one variable (env_var)
#> mutate (grouped): changed 672 values (100%) of 'Species' (0 new NA)
ggplot(wm, aes(year, `Weighted median`, color = Species, fill = Species)) +
geom_line() +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
ggplot(wm, aes(year, `Weighted median`, color = Species)) +
geom_line(size = 0.8) +
labs(linetype = "Quarter", x = "Year") +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
ggsave("figures/weighted_quantiles.pdf", width = 20, height = 20, units = "cm")
# Scaled
wm %>%
group_by(Species, quarter, env_var) %>%
mutate(`Weighted median scaled` = `Weighted median` - mean(`Weighted median`)) %>%
ggplot(aes(year, `Weighted median scaled`, color = Species)) +
geom_line(size = 0.8) +
labs(linetype = "Quarter", x = "Year") +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> group_by: 3 grouping variables (Species, quarter, env_var)
#> mutate (grouped): new variable 'Weighted median scaled' (double) with 647 unique values and 0% NA
ggsave("figures/supp/weighted_quantiles_scaled.pdf", width = 20, height = 20, units = "cm")
# ggplot(wm, aes(year, Median, color = Species, fill = Species)) +
# geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.4, color = NA) +
# facet_wrap(env_var ~ quarter, scales = "free", ncol = 2) +
# scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
# scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#
# ggsave("figures/supp/weighted_quantiles_range.pdf", width = 20, height = 20, units = "cm")
cpue3 <- cpue2 %>%
mutate(tot_cod = large_cod + small_cod) %>%
filter(quarter == 4)
#> mutate: new variable 'tot_cod' (double) with 7,870 unique values and 0% NA
#> filter: removed 5,271 rows (59%), 3,595 rows remaining
spde <- make_mesh(cpue3, xy_cols = c("X", "Y"),
n_knots = 200,
type = "kmeans", seed = 42)
cod_ind <- sdmTMB(tot_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = cpue3,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
fle_ind <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = cpue3,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
sanity(cod_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(fle_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# cod_ind2 <- sdmTMB(tot_cod ~ 0 + fyear + s(depth_sc),
# data = cpue3,
# mesh = spde,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# fle_ind2 <- sdmTMB(flounder ~ 0 + fyear + s(depth_sc),
# data = cpue3,
# mesh = spde,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# sanity(cod_ind2)
# sanity(fle_ind2)
Now try old data as well
# d_old <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
#
# # Filter to match the data I want to predict on
# d_old <- d_old %>%
# filter(year > 1992 & year < 2020 & quarter == 4)
#
# # Calculate standardized variables
# d_old <- d_old %>%
# mutate(depth_sc = (depth - mean(depth))/sd(depth),
# X = X/1000,
# Y = Y/1000,
# year = as.integer(year),
# fyear = as.factor(year)) %>%
# drop_na(depth) %>%
# rename("density_cod" = "density") # to fit better with how flounder is named
#
# nrow(d_old)
# nrow(cpue3)
#
# # Plot comparison
# p1 <- ggplot() +
# geom_histogram(data = d_old, aes(density_cod, fill = "old"), alpha = 0.5) +
# geom_histogram(data = cpue3, aes(tot_cod, fill = "new"), alpha = 0.5) +
# scale_y_continuous(trans = "log")
#
# p2 <- ggplot() +
# geom_histogram(data = d_old, aes(density_fle, fill = "old"), alpha = 0.5) +
# geom_histogram(data = cpue3, aes(flounder, fill = "new"), alpha = 0.5) +
# scale_y_continuous(trans = "log")
#
# p1 + p2
# spde_old <- make_mesh(d_old, xy_cols = c("X", "Y"),
# n_knots = 200,
# type = "kmeans", seed = 42)
#
# cod_ind3 <- sdmTMB(density_cod ~ 0 + fyear + s(depth_sc),
# data = d_old,
# mesh = spde_old,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# fle_ind3 <- sdmTMB(density_fle ~ 0 + fyear + s(depth_sc),
# data = d_old,
# mesh = spde_old,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# sanity(cod_ind3)
# sanity(fle_ind3)
### Cod
# SD 24
preds_cod24_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining
# SD 25
preds_cod25_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining
# SD 26
preds_cod26_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining
# SD 27
preds_cod27_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining
# SD 28
preds_cod28_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining
# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(4*4, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(4*4, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(4*4, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(4*4, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(4*4, nrow(preds_cod28_sim)))
# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
### Flounder
# SD 24
preds_fle24_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining
# SD 25
preds_fle25_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining
# SD 26
preds_fle26_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining
# SD 27
preds_fle27_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining
# SD 28
preds_fle28_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining
# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(4*4, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(4*4, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(4*4, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(4*4, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(4*4, nrow(preds_fle28_sim)))
# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
div_index_sim <- bind_rows(index24_cod_sim,
index25_cod_sim,
index26_cod_sim,
index27_cod_sim,
index28_cod_sim,
index24_fle_sim,
index25_fle_sim,
index26_fle_sim,
index27_fle_sim,
index28_fle_sim) %>%
mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes
#> mutate: new variable 'est_t' (double) with 280 unique values and 0% NA
#> new variable 'lwr_t' (double) with 280 unique values and 0% NA
#> new variable 'upr_t' (double) with 280 unique values and 0% NA
Plot the index and save as csv
# # First compare with previous index:
# cod_fle_index_old <- read.csv("/Users/maxlindmark/Dropbox/Max work/R/cod_condition/output/cod_fle_index.csv") %>%
# mutate(sub_div = as.character(sub_div)) %>%
# filter(sub_div %in% c("25", "28"))
#
# index_comp <- bind_rows(cod_fle_index_old %>% mutate(source = "cod_condition"),
# div_index_sim %>% mutate(source = "cod_interactions"))
#
# ggplot(index_comp, aes(year, est_t, color = source, fill = source)) +
# facet_wrap(sub_div ~ species, scales = "free") +
# geom_line() +
# geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t), alpha = 0.2, color = NA)
#
#
# ## If it isn't the same with the data, plot pred grids against each other
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2, color = NA) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
div_index_sim %>%
filter(sub_div %in% c(25, 28)) %>%
ggplot(aes(year, est_t/1000, color = species, fill = species)) +
geom_line() +
facet_wrap(~sub_div, scales = "free", ncol = 1) +
geom_ribbon(aes(year, ymin = lwr_t/1000, ymax = upr_t/1000), alpha = 0.2, color = NA) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "Species") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
#> filter: removed 168 rows (60%), 112 rows remaining
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
#geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
ggplot(div_index_sim, aes(year, est_t, color = species)) +
geom_line() +
facet_wrap(~sub_div, scales = "free") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
write.csv(div_index_sim, "output/cod_fle_index.csv")